## Authors: Alexander Herzog and Kenneth Benoit
## Date: May 30, 2015
## Replication file for JOP article "The Most Unkindest Cuts: Speaker Selection and Expressed Government Dissent During Economic Crisis"

rm(list = ls(all = TRUE))

library(rjags)
set.seed(42)

### PATHS ########################################
dataAll <- "./generated_data/working_data_all.RData"
##################################################

# MCMC values
mcmc.adapt <- 500
mcmc.burn <- 1000
mcmc.sample <- 2000
mcmc.thin <- 2
mcmc.chains <- 2

# JAGS model
JagsModel <- "model{
    for (i in 1:N) {
        
        # outcome model
        # -------------
        Yo[i] ~ dnorm(muo[i], tau)
        muo[i] <- ao[m[i]] + inprod(bo, Xo[i,1:Ko])

        # selection model, conditional on outcome model		
        # ---------------------------------------------
        Ys[i] ~ dbern(pis[i])
		
	# probit model for mapping normal latent variable on binary outcome		
        pis[i] <- phi(cmu[i]/csd[i])
        
	# calculate mean and precision, conditional on outcome model
	# see http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Bivariate_case
        cmu[i] <- mus[i] + (rho/sigma)*(Yo[i] - muo[i])
        csd[i] <- sqrt(1 - rho*rho)

        mus[i] <- as[m[i]] + inprod(bs, Xs[i,1:Ks])
    }
	
    # priors for selection coefficients
    # ---------------------------------
    for (j in 1:J)
    {
            as[j] ~ dnorm(mu.as, tau.as)
    }

    mu.as ~ dnorm(0, 0.0001)
    tau.as <- pow(sigma.as, -2)
    sigma.as ~ dunif(0, 1000)


    for (ks in 1:Ks) {
            bs[ks] ~ dnorm(0, 0.0001)
	}

    # priors for outcome coefficients	
    # -------------------------------
    for (j in 1:J)
    {
            ao[j] ~ dnorm(mu.ao, tau.ao)
    }

    mu.ao ~ dnorm(0, 0.0001)
    tau.ao <- pow(sigma.ao, -2)
    sigma.ao ~ dunif(0, 1000)

    for (ko in 1:Ko) {
        bo[ko] ~ dnorm(0, 0.0001)
    }	

    # priors for outcome variance and outcome-selection correlation
    # -------------------------------------------------------------
    tau ~ dgamma(0.01,0.01)
    sigma <- 1/sqrt(tau)
    rho ~ dunif(-1,1)
    
}"


# function for JAGS data for multilevel models
forJagsData <- function(DV.selection, DV.outcome, xMatrixSelection, xMatrixOutcome, groupID) {
    list(
        Ys = DV.selection,
        Yo = DV.outcome,
        Xs = xMatrixSelection,
        Xo = xMatrixOutcome,
        N = nrow(xMatrixSelection), 
        Ks = ncol(xMatrixSelection),
        Ko = ncol(xMatrixOutcome),
        J = length(unique(groupID)),
        m = groupID
    )
}

load(dataAll)


# Model 1 without interaction effects
# -----------------------------------
print("Estimating joint model 1")

selectionEq <- ~ 0 + lr_propScaled + safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + crisis
Xs <- model.matrix(selectionEq, dataAll)

outcomeEq <-  ~ 0 + lr_propScaled + safetyScaled + seniorityYearsScaled + government + cabMember + crisis
Xo <- model.matrix(outcomeEq, dataAll)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(dataAll$spoke, dataAll$textscore, Xs, Xo, dataAll$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

joint1.posterior <- coda.samples(model,
                                 variable.names=c("as", "ao", "bs", "bo", "sigma", "rho", "mu.as", "sigma.as", "mu.ao", "sigma.ao"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(joint1.posterior,file="./estimated_models/joint_model1_multilevel_posterior.RData")


# Model 2 with interaction effects
# ---------------------------------
print("Estimating joint model 2")

selectionEq <- ~ 0 + lr_propScaled + safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + crisis + crisis*lr_propScaled + crisis*safetyScaled
Xs <- model.matrix(selectionEq, dataAll)

outcomeEq <-  ~ 0 + lr_propScaled + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + crisis*lr_propScaled + crisis*safetyScaled
Xo <- model.matrix(outcomeEq, dataAll)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(dataAll$spoke, dataAll$textscore, Xs, Xo, dataAll$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

joint2.posterior <- coda.samples(model,
                                 variable.names=c("as", "ao", "bs", "bo", "sigma", "rho", "mu.as", "sigma.as", "mu.ao", "sigma.ao"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(joint2.posterior,file="./estimated_models/joint_model2_multilevel_posterior.RData")


# Model 3 with three-way interaction effects
# ------------------------------------------
print("Estimating joint model 3")

selectionEq <- ~ 0 + lr_propScaled + safetyScaled + seniorityYearsScaled + party.leader + government + cabMember + log.party.size + log.debate.days + crisis + crisis*lr_propScaled + crisis*safetyScaled + crisis*government + government*lr_propScaled + government*safetyScaled + crisis*government*lr_propScaled + crisis*government*safetyScaled
Xs <- model.matrix(selectionEq, dataAll)

outcomeEq <-  ~ 0 + lr_propScaled + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + crisis*lr_propScaled + crisis*safetyScaled + crisis*government + government*lr_propScaled + government*safetyScaled + crisis*government*lr_propScaled + crisis*government*safetyScaled
Xo <- model.matrix(outcomeEq, dataAll)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(dataAll$spoke, dataAll$textscore, Xs, Xo, dataAll$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

joint3.posterior <- coda.samples(model,
                                 variable.names=c("as", "ao", "bs", "bo", "sigma", "rho", "mu.as", "sigma.as", "mu.ao", "sigma.ao"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(joint3.posterior,file="./estimated_models/joint_model3_multilevel_posterior.RData")
